---
title: "APSR_COVID_code_LaFleur"
author: "LaFleur Stephens"
date: "2021/11/07"
output:
  pdf_document:
    keep_tex: true
header-includes:
   - \usepackage{dcolumn}    
classoption: landscape
fig_width: 6
fig_height: 4
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
rm(list=ls())
library(ggplot2)
library(stargazer)
library(sandwich)
library(dplyr)
library(pwr)
library(aod)
library(readstata13)
library(xtable)
library(arm)
library(psych)
library(tidyr)
library(purrr)
library(readr)
library(tibble)
library(stringr)
library(cobalt)
# Load data and create summary statistics
data <- read.dta13("/Users/Dropbox/Tess_Covid_update.dta")
data_pilot <- read.dta13("/Users/Dropbox/PilotData_June2021.dta")
```


# Table 1 - logit regression
```{r, results='asis'}
# Table 1: Influence of Racial Disparity Treatment and Negative Stereotype Endorsement
# Logit regressions for "negstereotype_endorsement"
glm1 <- glm(individualrights_dichotomous ~ Treatment + negstereotype_endorsement + Treatment * negstereotype_endorsement, family = binomial(link = "logit"), data = data, weights = WEIGHT)
glm2 <- glm(facemasks_dichotomous ~ Treatment + negstereotype_endorsement + Treatment * negstereotype_endorsement, family = binomial(link = "logit"), data = data, weights = WEIGHT)
glm3 <- glm(visitparks_dichotomous ~ Treatment + negstereotype_endorsement + Treatment * negstereotype_endorsement, family = binomial(link = "logit"), data = data, weights = WEIGHT)
glm4 <- glm(Follow_Blacks_dichotomous ~ Treatment + negstereotype_endorsement + Treatment * negstereotype_endorsement, family = binomial(link = "logit"), data = data, weights = WEIGHT)


stargazer(glm2, glm1, glm4, glm3, title = "Table 1: The Influence of Racial Disparities Treatment and Negative Stereotype Endorsement on COVID-19 Opinion", column.sep.width = "1pt", align = TRUE,  no.space=TRUE, dep.var.labels = c("Facemasks Are Not Important", "Individual Rights Threatened", "Black People Rarely Follow Social Distancing Guidelines", "Visit Parks Without Restrictions"), covariate.labels = c("Racial Disparities Information", "Negative Stereotype Endorsement", "Racial Disparities Information * \nNegative Stereotype Endorsement"), notes = c("Entries are logit coefficients. Standard errors in parentheses"), style = "apsr", out = "table1.pdf", type = "latex", model.numbers = F)
```

# Predicted probabilities - logit model
```{r, results='asis'}
# Below are code for predicted probabilities using logit model 
# Predicted probability "individualrights_dichotomous"
# Treatment group, negstereotype_endorsement = 1
p1.1 <- invlogit(coef(glm1)[1] + coef(glm1)[2] * 1 + coef(glm1)[3] * 1 + coef(glm1)[4] * 1)
# Treatment group, negstereotype_endorsement = 0
p2.1 <- invlogit(coef(glm1)[1] + coef(glm1)[2] * 1 + coef(glm1)[3] * 0 + coef(glm1)[4] * 0)
# Control group, negstereotype_endorsement = 1
p3.1 <- invlogit(coef(glm1)[1] + coef(glm1)[2] * 0 + coef(glm1)[3] * 1 + coef(glm1)[4] * 0)
# Control group, negstereotype_endorsement = 0
p4.1 <- invlogit(coef(glm1)[1] + coef(glm1)[2] * 0 + coef(glm1)[3] * 0 + coef(glm1)[4] * 0)

# Predicted standard errors
predict(glm1, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==1), se.fit=TRUE)
predict(glm1, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==1), se.fit=TRUE)
predict(glm1, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==0), se.fit=TRUE)
predict(glm1, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==0), se.fit=TRUE)
```


```{r, results='asis'}
# Predicted probability "facemasks_dichotomous"
  p1.2 <- invlogit(coef(glm2)[1] + coef(glm2)[2] * 1 + coef(glm2)[3] *1 + coef(glm2)[4])
  p2.2 <- invlogit(coef(glm2)[1] + coef(glm2)[2] * 1 + coef(glm2)[3] *0 + coef(glm2)[4] * 0)
  p3.2 <- invlogit(coef(glm2)[1] + coef(glm2)[2] * 0 + coef(glm2)[3] *1 + coef(glm2)[4] * 0)
  p4.2 <- invlogit(coef(glm2)[1] + coef(glm2)[2] * 0 + coef(glm2)[3] *0 + coef(glm2)[4] * 0)

predict(glm2, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==1), se.fit=TRUE)
predict(glm2, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==1), se.fit=TRUE)
predict(glm2, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==0), se.fit=TRUE)
predict(glm2, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==0), se.fit=TRUE)
```


```{r, results='asis'}
# Predicted probability "visitparks_dichotomous"
  p1.3 <- invlogit(coef(glm3)[1] + coef(glm3)[2] * 1 + coef(glm3)[3] *1 + coef(glm3)[4])
  p2.3 <- invlogit(coef(glm3)[1] + coef(glm3)[2] * 1 + coef(glm3)[3] *0 + coef(glm3)[4] * 0)
  p3.3 <- invlogit(coef(glm3)[1] + coef(glm3)[2] * 0 + coef(glm3)[3] *1 + coef(glm3)[4] * 0)
  p4.3 <- invlogit(coef(glm3)[1] + coef(glm3)[2] * 0 + coef(glm3)[3] *0 + coef(glm3)[4] * 0)

predict(glm3, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==1), se.fit=TRUE)
predict(glm3, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==1), se.fit=TRUE)
predict(glm3, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==0), se.fit=TRUE)
predict(glm3, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==0), se.fit=TRUE)

```


```{r, results='asis'}
# Predicted probability "followBlacks_dichotomous"
  p1.4 <- invlogit(coef(glm4)[1] + coef(glm4)[2] * 1 + coef(glm4)[3] *1 + coef(glm4)[4])
  p2.4 <- invlogit(coef(glm4)[1] + coef(glm4)[2] * 1 + coef(glm4)[3] *0 + coef(glm4)[4] * 0)
  p3.4 <- invlogit(coef(glm4)[1] + coef(glm4)[2] * 0 + coef(glm4)[3] *1 + coef(glm4)[4] * 0)
  p4.4 <- invlogit(coef(glm4)[1] + coef(glm4)[2] * 0 + coef(glm4)[3] *0 + coef(glm4)[4] * 0)

# Data frame with predicted probabilities
  allp4 <- rbind.data.frame(p1.4, p2.4, p3.4, p4.4)
  predict(glm4, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==1), se.fit=TRUE)
  predict(glm4, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==1), se.fit=TRUE)
  predict(glm4, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==0), se.fit=TRUE)
  predict(glm4, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==0), se.fit=TRUE)
```
  
  
# Figure 1 & Figure 2
```{r, results= 'asis'}
# Figure 1: Figure combining negstereotype=1 
all_lazy2_1 <- rbind.data.frame(p1.1, p3.1, p1.2, p3.2, p1.3, p3.3, p1.4, p3.4)

se_lazy2_1 <- rbind.data.frame(0.07744278, 0.07207574, 0.06431665, 0.0817368,0.06993488,0.04867954, 0.04160507, 0.06769423)
upper_lazy2_1 <- rbind.data.frame(se_lazy2_1*1.44 + all_lazy2_1)
lower_lazy2_1 <- rbind.data.frame(-se_lazy2_1*1.44 + all_lazy2_1)

all_condition <- rbind.data.frame("Treatment", "Control", "Treatment", "Control","Treatment", "Control", "Treatment", "Control")
  name_l <- rbind.data.frame("Individual Rights \nThreatened","Individual Rights \nThreatened", "Facemasks Are \nNot Important", "Facemasks Are \nNot Important", "Visit Parks \nWithout Restrictions", "Visit Parks \nWithout Restrictions", "Black People Rarely \nFollow Social \nDistancing Guidelines","Black People Rarely \nFollow Social \nDistancing Guidelines")
  
  all_lazy2_1 <- cbind.data.frame(name_l, all_lazy2_1, all_condition,upper_lazy2_1, lower_lazy2_1)

  colnames(all_lazy2_1) <- c("Variable", "Probability", "Treatment", "upper_1", "lower_1")
  
  ggplot(all_lazy2_1, aes(Treatment, Probability, fill = Treatment)) + geom_col() + ylab("Predicted Probability") + xlab("")+ theme(plot.title = element_text(hjust = 0.5, size = 11))+ theme(legend.title = element_blank()) +facet_wrap(~factor(Variable,levels = c('Facemasks Are \nNot Important', 'Individual Rights \nThreatened','Black People Rarely \nFollow Social \nDistancing Guidelines', 'Visit Parks \nWithout Restrictions')))+ scale_fill_manual("Group", values = c(Treatment = "grey36", Control = "grey"))  + geom_errorbar(aes(ymin=lower_1, ymax= upper_1), width=0.1)+ theme(plot.caption = element_text(hjust = 0)) 
  
  
  # Figure combining lazy2=0
  all_lazy2_0 <- rbind.data.frame(p2.1, p4.1, p2.2, p4.2, p2.3, p4.3, p2.4, p4.4)
  
  se_lazy2_0 <- rbind.data.frame (0.02615949, 0.03175135, 0.03104349, 0.03366404,0.02619647,0.0298043, 0.03199511,  0.03241079)
  upper_lazy2_0<- rbind.data.frame(all_lazy2_0 + 1.44* se_lazy2_0)
  lower_lazy2_0 <- rbind.data.frame(all_lazy2_0 - 1.44* se_lazy2_0)
  
  all_lazy2_0 <- cbind.data.frame(name_l, all_lazy2_0, all_condition, upper_lazy2_0, lower_lazy2_0)
    colnames(all_lazy2_0) <- c("Variable", "Probability", "Treatment", "upper_0", "lower_0")
    
  # Figure 2: Figure combining lazy2=0
 ggplot(all_lazy2_0, aes(Treatment, Probability, fill = Treatment)) + geom_col() + ylab("Predicted Probability") + xlab("") + theme(plot.title = element_text(hjust = 0.5, size = 11))+ theme(legend.title = element_blank()) +facet_wrap(~factor(Variable,levels = c('Facemasks Are \nNot Important', 'Individual Rights \nThreatened','Black People Rarely \nFollow Social \nDistancing Guidelines', 'Visit Parks \nWithout Restrictions')))+ scale_fill_manual("Group", values = c(Treatment = "grey36", Control = "grey")) + geom_errorbar(aes(ymin=lower_0, ymax= upper_0), width=0.1)+ theme(plot.caption = element_text(hjust = 0))
```


# Figure 3 -dif-in-dif_logit
```{r, results='asis'}
# Figure 3
### Difference in Difference with Dichotomized DVs
newdata1 = subset(data, negstereotype_endorsement==1 & Treatment==1)
newdata2 = subset(data, negstereotype_endorsement==0 & Treatment==1)
newdata3 = subset(data, negstereotype_endorsement==1 & Treatment==0)
newdata4 = subset(data, negstereotype_endorsement==0 & Treatment==0)

d1.1 <- (mean(newdata1$individualrights_dichotomous, na.rm= T)-mean(newdata3$individualrights_dichotomous, na.rm = T)) - (mean(newdata2$individualrights_dichotomous, na.rm = T)-mean(newdata4$individualrights_dichotomous, na.rm = T))
d2.1 <- (mean(newdata1$facemasks_dichotomous, na.rm= T)-mean(newdata3$facemasks_dichotomous, na.rm = T)) - (mean(newdata2$facemasks_dichotomous, na.rm = T)-mean(newdata4$facemasks_dichotomous, na.rm = T))
d3.1 <- (mean(newdata1$visitparks_dichotomous, na.rm= T)-mean(newdata3$visitparks_dichotomous, na.rm = T)) - (mean(newdata2$visitparks_dichotomous, na.rm = T)-mean(newdata4$visitparks_dichotomous, na.rm = T))
d4.1 <- (mean(newdata1$Follow_Blacks_dichotomous, na.rm= T)-mean(newdata3$Follow_Blacks_dichotomous, na.rm = T)) - (mean(newdata2$Follow_Blacks_dichotomous, na.rm = T)-mean(newdata4$Follow_Blacks_dichotomous, na.rm = T))

alld1 <- rbind.data.frame(d2.1, d1.1, d4.1, d3.1)

dlm1 <- lm(individualrights_dichotomous ~ Treatment + negstereotype_endorsement+ Treatment * negstereotype_endorsement, data = data)
dlm2 <- lm(facemasks_dichotomous ~ Treatment + negstereotype_endorsement+ Treatment * negstereotype_endorsement, data = data)
dlm3 <- lm(visitparks_dichotomous ~ Treatment + negstereotype_endorsement+ Treatment * negstereotype_endorsement, data = data)
dlm4 <- lm(Follow_Blacks_dichotomous ~ Treatment+ negstereotype_endorsement + Treatment * negstereotype_endorsement, data = data)

summary(dlm3)

upper_d <- c(d2.1 + 1.44 *0.12109, d1.1 + 1.44 * 0.11029, d4.1 + 1.44 *0.11618, d3.1 + 1.44 * 0.102517)
  
lower_d <- c(d2.1 - 1.44 *0.12109, d1.1 - 1.44 * 0.11029, d4.1 - 1.44 *0.11618, d3.1 - 1.44 * 0.102517)


named <- rbind.data.frame("Facemasks Are \nNot Important", "Individual Rights \nThreatened",  "Black People Rarely Follow \n Social Distancing Rules", "Visit Parks \nWithout Restrictions")

alld1 <- cbind.data.frame(named, alld1)
colnames(alld1) <- c("Group", "Probability")


ggplot(alld1, aes(Group, Probability)) +geom_col()  + ylab("")+ xlab("") + theme(plot.title = element_text(hjust = 0.5, size = 11))+ theme(legend.title = element_blank()) + theme(axis.text.x = element_text(size = 7.5))+ theme(plot.caption = element_text(hjust = 0)) +  scale_x_discrete(limits=c("Facemasks Are \nNot Important","Individual Rights \nThreatened",  "Black People Rarely Follow \n Social Distancing Rules", "Visit Parks \nWithout Restrictions")) + geom_errorbar(aes(ymin=lower_d, ymax= upper_d), width=0.1)



```


##### Code for appendix

# Table A2: pilot study
```{r, results='asis'}
# Regression with pilot data
lm1 <- glm(financiallypunished_dichotomous ~ newExperimentalConditions + newlazy_stereo+ newExperimentalConditions * newlazy_stereo, family = binomial(link = "logit"), data = data_pilot)

lm2 <- glm(stayathome_dichotomous ~ newExperimentalConditions + newlazy_stereo+ newExperimentalConditions * newlazy_stereo, family = binomial(link = "logit"), data = data_pilot)

lm3 <- glm(shutdown_dichotomous ~ newExperimentalConditions + newlazy_stereo+ newExperimentalConditions * newlazy_stereo, family = binomial(link = "logit"), data = data_pilot)

# Table A2
stargazer(lm1, lm2, lm3, title = "Table A2: The Influence of Racial Disparities Treatment and Negative Stereotype Endorsement on COVID-19 Opinion", covariate.labels = c("Racial Disparities Information", "Negative stereotype endorsement", "Racial Disparities Information * \nNegative Stereotype Endorsement"), dep.var.labels = c("Fining Risky Behavior", "Reopen the State \nProtestors Should Stay Home", "Non-Essential Business \nShould be Shut Down"), style = "apsr", notes = c("Entries are logit coefficients. Standard errors in parentheses"), model.numbers = F)
```
# Note that due to data sensitivity, no code nor data will be provided to reproduce Tables A3-A5


# Table A6: OLS Regressions
```{r, results='asis'}
# Regress non-dichotomized variables - OLS
m1 <- lm(individualrights ~ Treatment + negstereotype_endorsement + Treatment * negstereotype_endorsement, data = data, weights = WEIGHT)
m2 <- lm(facemasks ~ Treatment + negstereotype_endorsement + Treatment * negstereotype_endorsement, data = data, weights = WEIGHT)
m3 <- lm(Follow_Blacks ~ Treatment + negstereotype_endorsement + Treatment * negstereotype_endorsement, data = data, weights = WEIGHT)
m4 <- lm(visitparks ~ Treatment + negstereotype_endorsement + Treatment * negstereotype_endorsement, data = data, weights = WEIGHT)


stargazer(m2, m1, m3, m4, title = "Table A2: The Influence of Racial Disparities Treatment and Negative Stereotype \nEndorsement on COVID-19 Opinion (Non-dichotomized Dependent Variables)", column.sep.width = "1pt", align = TRUE,  no.space=TRUE, dep.var.labels = c("Facemasks Are Not Important", "Individual Rights Threatened", "Black People Rarely Follow Social Distancing Guidelines", "Visit Parks Without Restrictions"), covariate.labels = c("Racial Disparities Information", "Negative Stereotype Endorsement", "Racial Disparities Information * \nNegative Stereotype Endorsement"), notes = c("Entries are OLS coefficients. Standard errors in parentheses"), style = "apsr", out = "table1.pdf", type = "latex", model.numbers = F)
```

# Predicted Probability using OLS model
```{r, results='asis'}
# OLS individualrights predicted probability
summary(m1)
predict(m1, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==1), se.fit=TRUE)
predict(m1, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==1), se.fit=TRUE)
predict(m1, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==0), se.fit=TRUE)
predict(m1, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==0), se.fit=TRUE)

# OLS facemasks predicted probability
summary(m2)
predict(m2, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==1), se.fit=TRUE)
predict(m2, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==1), se.fit=TRUE)
predict(m2, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==0), se.fit=TRUE)
predict(m2, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==0), se.fit=TRUE)

# OLS Follow_Blacks predicted probability
summary(m3)
predict(m3, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==1), se.fit=TRUE)
predict(m3, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==1), se.fit=TRUE)
predict(m3, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==0), se.fit=TRUE)
predict(m3, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==0), se.fit=TRUE)

# OLS Visitparks predicted probability
summary(m4)
predict(m4, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==1), se.fit=TRUE)
predict(m4, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==1), se.fit=TRUE)
predict(m4, type="response", newdata = subset(data, negstereotype_endorsement==1 & Treatment==0), se.fit=TRUE)
predict(m4, type="response", newdata = subset(data, negstereotype_endorsement==0 & Treatment==0), se.fit=TRUE)


predicted_p_prejudiced <- rbind.data.frame(0.549, 0.487, 0.464, 0.364, 0.632, 0.544, 0.686, 0.557)
predicted_p_nprejudiced <- rbind.data.frame(0.355, 0.384, 0.3, 0.309, 0.587, 0.599, 0.444, 0.478)

se_p <- rbind.data.frame(0.0501189, 0.05211489, 0.04645526, 0.0487294, 0.03590808, 0.03945567, 0.03236762, 0.03404578)
se_np<- rbind.data.frame(0.02140435, 0.02254481, 0.02002452, 0.02113269, 0.01508053,0.01627657, 0.01401067, 0.01476208)

upper_ols_p <- predicted_p_prejudiced + 1.44 * se_p
lower_ols_p <- predicted_p_prejudiced - 1.44 * se_p
  
upper_ols_np<- predicted_p_nprejudiced + 1.44 * se_np
lower_ols_np<- predicted_p_nprejudiced - 1.44 * se_np


all_condition <- rbind.data.frame("Treatment", "Control", "Treatment", "Control", "Treatment", "Control","Treatment", "Control")
all_ols_p <- cbind.data.frame(name_l, predicted_p_prejudiced, all_condition, upper_ols_p, lower_ols_p)

all_ols_np <- cbind.data.frame(name_l, predicted_p_nprejudiced, all_condition, upper_ols_np, lower_ols_np)

colnames(all_ols_p) <- c("Variable", "Probability", "Treatment", "upper_ols_p", "lower_ols_p")
colnames(all_ols_np) <- c("Variable", "Probability", "Treatment", "upper_ols_np", "lower_ols_np")

# Predicted Probability with prejudiced
ggplot(all_ols_p, aes(Treatment, Probability, fill = Treatment)) + geom_col() + ylab("Predicted Probability") + xlab("")+ theme(plot.title = element_text(hjust = 0.5, size = 11))+ theme(legend.title = element_blank()) +facet_wrap(~factor(Variable,levels = c('Facemasks Are \nNot Important', 'Individual Rights \nThreatened','Black People Rarely \nFollow Social \nDistancing Guidelines', 'Visit Parks \nWithout Restrictions')))+ scale_fill_manual("Group", values = c(Treatment = "grey36", Control = "grey"))  + geom_errorbar(aes(ymin=lower_ols_p, ymax= upper_ols_p), width=0.1)+ theme(plot.caption = element_text(hjust = 0))

# Predicted Probability with non-prejudiced
ggplot(all_ols_np, aes(Treatment, Probability, fill = Treatment)) + geom_col() + ylab("Predicted Probability") + xlab("")+ theme(plot.title = element_text(hjust = 0.5, size = 11))+ theme(legend.title = element_blank()) +facet_wrap(~factor(Variable,levels = c('Facemasks Are \nNot Important', 'Individual Rights \nThreatened','Black People Rarely \nFollow Social \nDistancing Guidelines', 'Visit Parks \nWithout Restrictions')))+ scale_fill_manual("Group", values = c(Treatment = "grey36", Control = "grey"))  + geom_errorbar(aes(ymin=lower_ols_np, ymax= upper_ols_np), width=0.1)+ theme(plot.caption = element_text(hjust = 0))


# Figure for Non-dichotomized Difference in Difference (OLS)

d1.2 <- (mean(newdata1$individualrights, na.rm= T)-mean(newdata3$individualrights, na.rm = T)) - (mean(newdata2$individualrights, na.rm = T)-mean(newdata4$individualrights, na.rm = T))
d2.2 <- (mean(newdata1$facemasks, na.rm= T)-mean(newdata3$facemasks, na.rm = T)) - (mean(newdata2$facemasks, na.rm = T)-mean(newdata4$facemasks, na.rm = T))
d3.2 <- (mean(newdata1$visitparks, na.rm= T)-mean(newdata3$visitparks, na.rm = T)) - (mean(newdata2$visitparks, na.rm = T)-mean(newdata4$visitparks, na.rm = T))
d4.2 <- (mean(newdata1$Follow_Blacks, na.rm= T)-mean(newdata3$Follow_Blacks, na.rm = T)) - (mean(newdata2$Follow_Blacks, na.rm = T)-mean(newdata4$Follow_Blacks, na.rm = T))

alld2 <- rbind.data.frame(d2.2, d1.2, d4.2, d3.2)

summary(m4)

upper_d2 <- c(d2.2 + 1.44 *0.076527, d1.2 + 1.44 * 0.08208, d4.2 + 1.44 *0.05340, d3.2 + 1.44 * 0.06032)
lower_d2 <- c(d2.2 - 1.44 *0.076527, d1.2 - 1.44 * 0.08208, d4.2 - 1.44 *0.05340, d3.2 - 1.44 * 0.06032)
alld2 <- cbind.data.frame(named, alld2)
colnames(alld2) <- c("Group", "Probability")


ggplot(alld2, aes(Group, Probability)) +geom_col() + ylab("")+ xlab("") + theme(plot.title = element_text(hjust = 0.5, size = 11))+ theme(legend.title = element_blank()) + theme(axis.text.x = element_text(size = 7.5))+ theme(plot.caption = element_text(hjust = 0)) +  scale_x_discrete(limits=c("Facemasks Are \nNot Important","Individual Rights \nThreatened",  "Black People Rarely Follow \n Social Distancing Rules", "Visit Parks \nWithout Restrictions")) + geom_errorbar(aes(ymin=lower_d2, ymax= upper_d2), width=0.1)
```


#Figures A2-A5 # Distribution of four dependent variables
```{r, results='asis'}
# Distribution of Non-dichotomized dependent variables
ggplot(data = data, aes(facemasks)) + geom_bar()
ggplot(data = data, aes(individualrights)) + geom_bar()
ggplot(data = data, aes(Follow_Blacks)) + geom_bar()
ggplot(data = data, aes(visitparks)) + geom_bar()
```
#Figure A1: # Balance tests
```{r, results='asis'}
# Balance tests
t.test(AGE ~ Control, data)
d_age = 50.84488-51.09375
t.test(gender_new ~ Control, data)
d_gender = 0.5247525 - 0.4930556
t.test(ideo5 ~ Control, data)
d_ideo5 = 0.4653285 - 0.4637405
t.test(pid7 ~ Control, data)
d_pid7 =  0.5280641 - 0.5218978
t.test(educ ~ Control, data)
d_educ = 0.7271727 - 0.7233796
t.test(statemention ~ Control, data)
d_statemention = 0.1485149 - 0.1284722
t.test(South ~ Control, data)
d_South = 0.3003300- 0.3159722
t.test(INCOME ~ Control, data)
d_income = 10.60726-10.39236

dif_covariate <- rbind.data.frame(d_age, d_gender, d_ideo5, d_pid7, d_educ, d_statemention, d_South, d_income)
name_covairate <- rbind.data.frame("Age", "Gender", "Ideology", "Party Affiliation", "Education", "State Mentioned", "South", "Income")
lower_covariate <- rbind.data.frame(-3.063170, -0.04920263, -0.04676515, -0.04701312, -0.03802174, -0.03581462, -0.09040084, -0.4272374)
upper_covariate <- rbind.data.frame(2.565439, 0.11259647, 0.04994117, 0.05934578, 0.04560792, 0.07589988, 0.05911646, 0.8570366)
covariate <- cbind.data.frame(dif_covariate, name_covairate, lower_covariate, upper_covariate)
colnames(covariate) <- c("Difference", "Variable", "lower95", "upper95")
ggplot(data = covariate, aes(x= Variable, y= Difference))+ geom_point()+ geom_errorbar(aes(ymin=lower95, ymax=upper95),  width=0.2)+ ylab("Difference in Means")
```

From the above results, I would consider the covariates to be balanced across control and treatment groups, given that the means are similar and the p-value is not small enough to reject the null hypothesis (true different in means is equal to 0). By default, the t-test used here is Welch t-test (without assuming equal variances)
